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A model of clustering dynamics is proposed for a population of spatially distributed active rotators. 
A transition from excitable to oscillatory dynamics is induced by the increase of the local density of 
active rotators. It is interpreted as dynamical quorum sensing. In the oscillation regime, phase waves 
propagate without decay, which generates an effectively long-range interaction in the clustering 
dynamics. The clustering process becomes facilitated and only one dominant cluster appears rapidly 
as a result of the dynamical quorum sensing. An exact localized solution is found to a simplified 
model equation, and the competitive dynamics between two localized states is studied numerically. 
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Synchronization in populations of coupled oscillators have been extensively studied [T!] since the pioneering works 
by Winfree pj] and Kuramoto fs*]. In spatially extended systems such as the BZ reaction, target patterns and spiral 
patterns appear as a result of synchronization On the other hand, populations of unicellular organisms sometimes 
undergo drastic transitions when the cell density increases. For example, some kinds of bacteria change to produce 
toxic substances when the cell density is beyond a threshold. The phenomenon is called quorum sensing in the 
meaning that the transition occurs by sensing the density of the same kind of organisms. Dynamical quorum sensing 
is used for a phenomenon that populations of active elements undergo some dynamical transitions with increasing the 
density. Dynamical quorum sensing was studied for synthetic multi-cellular clock Q , cell-density-dependent glycolytic 
oscillation in yeast Q, and chemical oscillators Gregor et al. studied collective behaviors of social amoebae: 
Dictyostelium discoideum, and found a transition from an excitable state to an oscillatory one with increasing the 
cell density[8j. The oscillation of extra cellular cAMP is observed in the experiment. The social amoebae begin 
to aggregate when the period of the oscillation becomes shorter. Some local aggregates are created initially, but a 
competition among local aggregates occurs via cell-cell signaling. The cell-cell signaling is caused by the waves of 
cAMP. Finally, only one dominant aggregate appears. 

There is a mathematical model for the oscillation of cAMP by Martiel and Goldbeter P. However, Gregor et al. 
used a simpler phase equation to explain the transition from an excitable state to an oscillatory state. We use a 
similar type of phase equation which can express a transition from excitatory dynamics to oscillatory one: 

— = a; — 6 sin (A. (1) 
dt ^ ' 

When uj < h, there is a stable stationary solution (p = sin^^(cj/6). The stationary state can be regarded as excitable 
because it is easily excited by a small external input when w is close to b. When w > 6, the phase increases 
monotonically, which represents the limit-cycle oscillation. Gregor et al. studied a globally coupled system of the 
form: 



dt 



V K + CEf=i{(-A,„ + A) sin,/), + A„ + ; ■ 



In this paper, we consider a spatially-extended oscillatory medium in one dimension. We use a phase equation of the 

form 



-=a.-6sin,/>-f.— , (3) 

because this type of phase equation is more general, although Eq. (2) or its generalization might be better for the 
quantitative argument. The third term expresses the phase diffusion, and the last nonlinear term g{d4>/ dx)^ plays 
an important role in the wave propagation such as target and spiral patterns especially in oscillatory systems. The 
derivation of the diffusion term and the nonlinear term for general oscillatory systems was discussed in Ref . [3] . We use 
a term "active rotator" in this paper which exhibits excitable or oscillatory dynamics and moves actively by sensing 
its density. We assume that u) in Eq. (3) is a function of n as w = an/ (1 -I- ^n), taking into account an experimental 
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fact that the frequency increases with the cell density of the social amoebae and is saturated when the density is 
sufficiently large. 

In this paper, we would like to study interaction between the dynamical quorum sensing and the agg regation 
dynamics. The Keller-Segel model is a model of the aggregation dynamics induced by chemotaxis 1 111 ]. In the 
chemotaxis, the organisms move, sensing the concentration of chemical substances. In a variant of the Keller-Segel 
model, the density n{x,t) obeys 



dt dx^ dx \ dx 

where D is the diffusion constant for n and represents a sensitivity function for the concentration v of the chemical 
substance. The organism moves toward the region of high x(^)- Oscillatory dynamics is not explicitly taken into 
consideration in the Keller-Segel model. Van Oss et al. studied a complicated model combining the Martiel-Goldbcter 
model of the cAMP oscillation with a certain aggregation dynamics, and succeeded in reproducing the motion toward 
the center of a spiral pattern, i.e., the source of the cAMP waves. 

In this paper, we use a variant of Eq. (4) for the aggregation dynamics. In our model, x(^) Eq. (4) is replaced 
by (j). It is because the cells move in the opposite direction to the propagation of cAMP waves. As a result they tend 
to aggregate toward the source of cAMP waves. The waves are sent out from a pacemaker region where the phase ip 
is in advance. 

The purpose of this paper is to show that the aggregation dynamics changes qualitatively by the dynamical quorum 
sensing. Our model equations are expressed as 

d(j} an , ■ , , d'^'t' , f 

^sin(/' + t^Tr:7+5 TT- , (5) 



dt 1 + 7n dx^ \ dx 

dn ^d'^n d f dip \ 

DtH^ - it- ■ (6) 



dt dx"^ dx \dx 

We study the coupled model equations numerically and theoretically. For numerical simulations, periodic boundary 
conditions 0(0) = 4>{L) and n(0) = n{L) are imposed. The total number of active rotators N = n{x,t)dx 
is conserved in the time evolution of Eq. We will study the model equations (5) and (6) to understand the 
interaction between the dynamical quorum sensing and the aggregation dynamics qualitatively in this paper. However, 
it is noted that the model equations (5) and (6) have a limit of application for quantitative argument: the amplitude 
of the oscillation is assumed not to change vastly in the phase description, but the oscillatory dynamics of cAMP 
waves may change qualitatively when the density n{x,t) becomes quite low, and the phase description might become 
worse. 

There is a uniform and stationary solution 4'{x, t) = 0o and n{x, t) — uq ^ N/L, which satisfies 

— = 6sm0o, (7) 

1 + 7?io 

if 6 > ano/{l + jng). Small perturbations 54> — 4> — 4>q and 6n — n{x, t) — uq obey linear equations: 

aSn , , ^ , d^S4> 
ocos (pctoq) + v- 



dt (1 + 7no)^ dx'^ 
dSn ^d'^Sn d / dSip \ 
dt dx"^ dx \ dx °/ 



By the Fourier transform, the Fourier amplitudes Scjjk and duk with wavenumber k obey 

d6(l)k aSuk 



dt (1 + jno) 

dSuk 
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bcos(f>oS(f>k — vk Stpk, 



dt 

The eigenvalue Xk is evaluated at 



Dk'^dnk +nok^(j)k. (9) 
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-{Dk^ +vk'^ + b cos M + JiDk'^ - i/fc2 _ b cos cjyo)'^ + Aanok'^ / {I + juq) ,4 
Xk = oifc +a2k, (10) 
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FIG. 1: (a) Eigenvalue X{k) for a = 0.15. Time evolutions of n{x,t) at (b) a = 0.15, (c) 0.16, and (d) 0.18 for 7 = 0.5,!^ = 
0.1, g = 1, & = 0.2, D = 0.1 and L = 60. 



where 

^ ariQ a^rig (D — v)anQ 

'^^ 5cos0o(l + 7710)^ ' ^ (1 + 77io)''6'^ cos^ (/)o (1 + 7no)^6^ cos^ (/)o 

for sufficiently sniall k. The uniform state becomes unstable for D < ano/{6cos0o(l + T'^o)^}? a-nd active rotators 
begin to aggregate locally. 

Figure 1(a) shows the eigenvalue A(fc) by Eq. (10) at a = 0.15. The other parameters are fixed to be 7 = 0.5, v = 
0.1, g — 1, b = 0.2, D — 0.1, no = 1 and L — 60. The eigenvalue takes a maximum near k — km = 1.25. Figures 
l(b)-(d) show time evolutions of n{x,t) at (b) a = 0.15 and (c) a = 0.16 and (d) a ~ 0.18. The initial conditions are 
(j){x, 0) = and n(x, 0) = 1 + r{x) where r is a random number between -0.02 and 0.02. The uniform state is unstable, 
a wavy pattern appears. The number of peaks near t — 100 at a = 0.15 is 12, and therefore the average wavelength is 
evaluated at 60/12 — 5, which is close to 27r/fcm — 5.02 for a — 0.15 shown in Fig. 1(a). The linear stability analysis 
explains the length scale of initial wavy patterns. The wavy pattern develops into localized clusters of active rotators. 
The clusters are mutually merged in 100 < t < 200, and the cluster number decreases during the time. The localized 
clusters become almost stationary at a = 0.15. However, a rapid merging process among clusters occurs at i > 300 
for a = 0.16. We have shown only for t < 500 in Fig. 1(c), however, the merging process continues and one cluster 
appears finally. The merging process occurs more strongly at a = 0.18, and only one cluster dominates at i = 500. 
Figure 2(a) shows time evolutions of 4>{xo) at xq = 15 for a = 0.15, 0.16 and 0.18. The phase (j){xo) remains to be 
below tt/2 for a fairly long initial time, but it begins to increase abruptly a.t t ^ 300 for a = 0.16 and at i ~ 140 
for a — 0.18, which implies that a transition from an excitable state into an oscillatory state has set in at the times. 
Stepwise time evolutions are seen for a = 0.16 and 0.18 in the oscillatory state. It is because the time derivative 
d(f>/dt becomes small when <j> goes through 7r/2 + 2?™. Figure 2(b) shows the time evolution of the profile (f>{x,t) at 
a = 0.16. The overlapped profiles below <j) — 7r/2 correspond to 0(x)'s for t < 300. The phase patterns grow abruptly 
at t ~ 300 by the transition. The phase profile has a peak around x ^ xq — 15, which implies that phase waves are 
sent out from a; ~ xq. Because the active rotators move to the center of the phase waves in our model, the localized 
clusters move toward the center and they are absorbed into a dominant cluster near x = xq at a = 0.16 as shown 
in Fig. 1(c). At a = 0.15, the transition to the oscillatory state does not occur and therefore the rapid merging of 
localized clusters is not observed. We can interpret the transition from the excitable state to the oscillatory state as 
a kind of dynamical quorum sensing in that the transition occurs when the local density increases and goes beyond a 
threshold. The clustering process changes qualitatively by the dynamical quorum sensing. In the excitatory regime, 
the clustering proceeds locally, however, a rapid clustering into one dominant cluster does not occur. This type of 
local clustering is similar to that in the Keller-Segel models [l^, in that the oscillatory dynamics is not involved. In 
the oscillation regime, localized clusters are attracted to the center of phase waves and are absorbed into the dominant 
cluster even if the localized clusters are far distant from the dominant cluster. It is because the phase waves propagate 
far away without decay in the oscillation regime. 

The transition from the excitable state to the oscillatory state can be studied more easily for a one-hump state. In 
the excitatory regime, there is a localized stationary solution, which has a form: 

0(a;) — - ax^, n{x) — UsC''^^^ , (11) 
where 0s, a, Ug and c are parameters to characterize the form of the one-hump solution. Substitution of Eq. (|lip into 
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FIG. 2: (a) Time evolutions of (j}{x,t) at a:: = xq = 15 for a = 0.15,0.16 and 0.18. (b) Time evolution of 4i{x,t) at a = 0.16. 
The other parameters are 7 = 0.5, ly = 0.1, g — 1, b — 0.2, D = 0.1 and L = 60. 
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FIG. 3: (a) Stationary solutions n{x,t) (solid curve) by direct numerical simulation of Eqs. (5) and (6), and n^e"'^^ (dashed 
curve) for a = 0.16, ^ = 1, g — l,u — 0.1, D = 0.1 and N = 7. (b) Relations of the peak value (ps and A^. Solid curve shows 
numerical results by Ea. (ll2p . and rhombi denote direct numerical results for a — 0.16,7 — 0.5, g — 1,1/ = 0.1 and D — 0.1. 
The dashed curve denotes a relation of the maximum value (j)m of 4>{x,t) in 13 < t < 17 and A'^; = J^^ n{x,t)dx in the time 
evolution of direct numerical simulation shown in Fig. 1(c). 



Eq. (5) yields 

6sin0s — 2;/a = 0, ^ — — 6a cos — 4ga =0, (12) 

1+7^3 (l + 7?is)"^ 

by comparing the coefficients of the first and second terms in the Taylor expansion by x around x = 0. Furthermore, 
Dc = a is satisfied from the relation Ddn/dx = {d<j)/dx)n in the stationary state. The total number is expressed 
as = J^^n{x)dx — Us^/ttJc. Then, the parameters a and c are expressed as a = imlD/N'^ and c = T:nl/N^. 
Substitution of these relations into Eq. (fT2)) yields coupled equations for ips amd . We can solve the coupled equations 
numerically. The solution is Ug = 4.92 and c — 1.55 for a = 0.16,7 = l,g = l,v = 0.1, D = 0.1 and N — 7. Figure 
3(a) shows n{x,t) (soHd curve) in the stationary state by direct numerical simulation of Eqs. (5) and (6) and rt^e"'^^ 
(dashed curve) with Us = 4.92 and c = 1.55 for the same parameters: a — 0.16,7 = l,.g = l,t^ = 0.1,1? = 0.1 and 
N = 7. Fairly good agreement is seen. Figure 3(b) shows a relation of 4>s and N for a — 0.16,7 — 0.5, 5 = \,v — 0.1 
and _D = 0.1 by the direct numerical simulation (rhombi) and by numerical results (solid curve) using Eq. (jl2p . The 
two values of 4>s{N) are in good agreement. The stationary solution (ps increases with N and reaches 7r/2, and then 
the stationary solution disappears when N is greater than the critical value. It implies the transition to the oscillatory 
state. 

In the time evolution shown in Fig. 1(c), there is a localized structure of n{x,t) around x — xq ^ 15. The local 
density n{x, t) increases in time there. We have calculated the number of active rotators in the localized structure as 
Ni ~ J^^ n{x,t)dx where xi — 13 and X2 = 17 and the maximum value 0m of 4>(x,t) in 13 < x < 17. Ni increases 
monotonically, and 0m also increases with Ni. The two values reach Ni = 9.05 and (pm = 1-57 dX t = 267, and then 
the transition into an oscillatory state occurs. A relation of Ni and the phase 0m in the time evolution is shown 
in Fig. 3(b) by a dashed curve. The trajectory of {Ni,(p„i) follows the curve of the stationary solution for < 6, 
which implies the localized structure is close to the stationary solution, however, the trajectory deviates from that of 
one-hump solutions for > 6, and iV; increases to 9.05 and then a transition to the oscillatory state occurs. The 
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FIG. 4: (a) Time evolution of n{x,t) for 6 = 0, 7 = 0, a = 0.5, g — 1 and A'' = 12. (b) Stationary solutions obtained by the 
direct numerical simulation (solid curve) and theoretical one (dashed curve). The difference is invisible. 
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FIG. 5: (a) Time evolution of cj){x,t). The initial conditions are n{x,0) — Ai/ cosh^ {ki{x — xi)} + A2/ cosh^{fc2(a; — 0:2)} where 
Ai = 2.0833, fei — 0.8333,2:1 — 18 and A2 = 1.333, fc2 = 0.666, and 12 = 42 (b) Time evolution of n{x,t) in the region of 
10 < a; < 50 for 50 < t < 90. (c) Time evolution of the peak position Xp2 of the second localized state. The dashed line is a 
line of velocity —0.833. (d) Time evolution of the peak position x-p2 of the second localized state for b = 0.25 and 7 = 0. The 
dashed line is a line of velocity —0.45. 



deviation might be due to the non-stationarity or the rapid increase of Ni. In any case, the critical value of 4>m is 
around (pm ~ 7r/2 = 1.57 even for this non-stationary localized structure. 

In the oscillation regime, many localized structures compete with each other, and only one localized structure 
survives finally. The dynamics in the oscillation regime can be studied using a simpler system of & = and 7 = 0. 
Figure 4(a) shovifs the time evolution of n{x, t) for a = 0.5, v = 0.5, g = 1,D = 0.5 and N = 12. The initial condition 
is (^(a:, 0) = and n{x,t) = 0.2 + r{x) where r is a random number between -0.001 and O.OOI. Three localized 
structures are created initially, and only one cluster survives by the competition among the three localized structures. 
Figure 4(b) shows the stationary solution localized around x = 22.3. The stationary solution can be explicitly solved 
in case of 6 = and 7 = by the ansatz n(x) = A/ cosh^(fca;) and (t>{x) — 13 In n + cot. Substitution of the ansatz into 
Eqs. (5) and (6) determines the parameter A, k and uj as 

4iyL' + 8(7132 ' SiyD + lGgD^' {2vD + AgDf ^ ' 

This solution is an exact solution of the nonlinear equations (5) and (6). The parameters are determined as fc = 2 
and A = 12 for a = 0.5,(7 — 1 ^nd N = 12. The dashed curve in Fig. 4(b) is the theoretical one: n{x) = 
12/cosh^{2(a; - 22.3)2}. The difference is invisible. 

To understand the competitive dynamics between two clusters, we have performed a numerical simulation from an 
initial conditions: (^(x,0) = and ri(x,0) = Ai/cosh^ ki{x — xi) + A2/ cosh^ k2{x — X2) where Ai = 2.0833, fci — 
0.8333, xi = 18 and A2 = 1.333,^2 = 0.666, 2:2 — 42. The parameters Ai,ki,A2 and ^2 are theoretical values by 
Eq. (12) respectively for A^i = 5 and iV2 = 4 at a = 0.5, v = D = 0.5 and g = 1- Figure 5(a) shows time evolution of 
(j){x,t). Initially, phase waves are sent out from the two center xi = 18 and X2 = 42. The phase waves from the two 
centers are expressed as 



(pi = D\nni{x) + LUit = D \n{Ai / cos'Ja (hx)} + Uit ±2DkiX + Dln{4:At) + ujit, for x Too, (14) 
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for i = 1 and 2. The two phase waves are approximated at ^ —2Dki{x — xi) + cjit and (f)2 ^ 2Dk2(x — x-i) + uj^t 
in the region of < x < x^- The region x\ < x < Xi separated into two regions by a boundary or a shock where 
the two phase waves collide with each other. The two localized structures arc fairly stable when the boundary is far 
from the two centers. However, the frequency Wi of the left waves is faster and therefore the boundary (the shock) 
between the two regions moves in the right direction. The boundary point is evaluated as x = (wi — ui2)t/{2D{ki + 
fe)} + (^22:2 + kixi)/{ki + k2) from the relation (f)i = (j)2, and therefore the velocity of the boundary is evaluated as 
V = (wi LJ2) /{2£'(fci + fe)}- When the boundary reaches X2, the localized structure of the right cluster begins to be 
destroyed. Figure 5(b) shows the time evolution of n{x, t) in the region of 10 < x < 50 for 50 < t < 90. The peak of 
the second localized structure moves to the left, the peak height decreases, and finally the right cluster is swallowed 
by the left dominant cluster. Figure 5(c) shows the peak position of the second (right) localized structure. The peak 
position is almost stationary until t = 70, but moves to the left for t > 72. The dashed line shows x ~ — 0.833f, 
that is, the velocity of the second peak is approximately —2Dki ~ —0.833. It is because the time evolution of n{x, t) 
around the second peak is approximated by 

dn j-)^^^ d4>dn d'^cf) d'^^^ +2Dk (15) 
dt dx"^ dx dx dx^ dx"^ dx ' 

owing to d(p/dx ^ —2Dki. Here, the phase gradient d(f)/dx is approximated at the phase gradient —2Dki of the phase 
wave sent out from the left center as is seen from Fig. 5(a). The last term 2Dki{dn/ dx) induces the motion of the 
second peak with velocity —2Dk\, and the first term makes the localized structure diffuse. Finally, a localized state 
appears around x = x\ with total number A/' = 4 + 5 = 9. Phase waves with a constant phase gradient propagate to 
the regions far away from the center, which is an origin of the effective long-range interaction in the oscillation regime. 
The effective long-range interaction in the oscillation regime facilitates the clustering process, and the clusters which 
are initially distant from one another are merged into one cluster rapidly. 

The approximation of the average velocity at —2Dki of the second peak is valid even for nonzero 6 if 6 is small, 
however, this approximation becomes worse near the transition to the excitable system. Figure 5(d) shows the peak 
(maximum) position of the second localized structure at b = 0.25 for 7 = 0, A''i = 5, 7V2 = 4, a = 0.5, v = D = 0.5 
and g = 1. The second localized structure is destroyed, a small two-peak structure appears, the two-peak structure 
exhibits oscillation, and moves toward the first localized structure on the average. The average velocity V2 of the 
two-peak structure is 0.54, however, the average phase gradient ki is about 0.9, and the equality V2 = —2Dk\ is not 
satisfied at this parameter. The discontinuous time evolution in Fig. 5(d) is due to the two-peak structure. 

In summary, we have proposed a model of clustering dynamics of spatially-distributed active rotators. A transition 
from an excitable state to an oscillatory state is induced by the increase of local density, which is interpreted as 
the dynamical quorum sensing. The transition occurs when the peak value of 4'{x,t) goes over 7r/2. The clustering 
dynamics changes qualitatively by the dynamical quorum sensing. That is, the clustering process becomes rapid, and 
only one cluster survives after the competition among localized clusters. Wc have found an exact localized sohition 
for a simpler system of 7 = 6 = 0, and studied the competitive dynamics between two localized states. In this paper, 
we have used a rather simplified model to understand the dynamics qualitatively. We would like to study a more 
quantitative model equation in the future, using a modified model of Eq. (2) and taking into account experimental 
results. It is another problem to extend the one-dimensional model equations (5) and (6) to a two-dimensional system. 
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